1 Summary

1.1 Groups information

mae <- get(load("hg38_gcimp_mae.rda"))
mae

## A MultiAssayExperiment object of 2 listed
##  experiments with user-defined names and respective classes.
##  Containing an ExperimentList class object of length 2:
##  [1] DNA methylation: RangedSummarizedExperiment with 135374 rows and 244 columns
##  [2] Gene expression: RangedSummarizedExperiment with 57035 rows and 244 columns
## Features:
##  experiments() - obtain the ExperimentList instance
##  colData() - the primary/phenotype DataFrame
##  sampleMap() - the sample availability DataFrame
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment
##  *Format() - convert into a long or wide DataFrame
##  assays() - convert ExperimentList to a SimpleList of matrices

samples.info <- SummarizedExperiment::colData(mae)
plyr::count(samples.info$subtype_Supervised.DNA.Methylation.Cluster)

1.2 Samples

1.3 Arguments

1.4 Summary results

1.4.1 Probes hypermethylated in G-CIMP-high vs G-CIMP-low

2 Glioma analysis: G-CIMP-high vs G-CIMP-low

2.1 Probes hypermethylated in G-CIMP-high vs G-CIMP-low

Volcano plot - Probes hypermethylated in G-CIMP-high vs G-CIMP-low

2.1.1 Significant anti-correlated pairs of gene-probes

Heatmap: hypomethylated paired probes

2.1.1.1 Statehub: Chromatin state evaluation

The plot below was produced with funciVar tool (https://github.com/Simon-Coetzee/funcivar), which calculate overlaps and enrichment between genomic variants and genomic features or segmentations. The segmentations used were retrieved from (http://statehub.org/).

Funcivar: hypermethylated paired probes. Produced with funciVar tool (https://github.com/Simon-Coetzee/funcivar) and statehub.org data.

2.1.2 Motif enrichment analysis

Probes hypermethylated in G-CIMP-high vs G-CIMP-low

2.1.3 TF analysis

2.1.3.1 Top potential TF

The table below shows for a given enriched motif, the top potential TF (best ranked TF based on the p-value) belonging to the same family or subfamily (TFClass classification) of the TF motif. The columns with a prefix p-value shows how significant is the anti-correlation of the average DNA methylation level of probes with the given motif and the TF expression.

2.1.3.2 potential TF family

The table below shows for a given enriched motif, all potential TF belonging to the same family (TFClass classification) of the TF motif. The columns with a prefix p-value shows how significant is the anti-correlation of the average DNA methylation level of probes with the given motif and the TF expression.

2.1.3.3 potential TF subfamily

The table below shows for a given enriched motif, all potential TF belonging to the same subfamily (TFClass classification) of the TF motif. The columns with a prefix p-value shows how significant is the anti-correlation of the average DNA methylation level of probes with the given motif and the TF expression.

2.1.3.4 TF plots

2.2 Probes hypomethylated in G-CIMP-high vs G-CIMP-low

Volcano plot - Probes hypomethylated in G-CIMP-high vs G-CIMP-low

3 Code

library(stringr)
library(TCGAbiolinks)
library(dplyr)
library(ELMER)
library(parallel)
#-----------------------------------
# 1 - Identify Samples
# ----------------------------------
gcimp.samples <- TCGAquery_subtype("lgg.gbm") %>% filter(base::grepl("G-CIMP",Supervised.DNA.Methylation.Cluster,ignore.case = T))
gcimp.low.samples <- gcimp.samples %>% filter(base::grepl("G-CIMP-low",Supervised.DNA.Methylation.Cluster, ignore.case = T))
gcimp.high.samples <- gcimp.samples %>% filter(base::grepl("G-CIMP-high",Supervised.DNA.Methylation.Cluster, ignore.case = T))

msg <- paste0("-----------------------------------------------\n",
              "GCIMP samples ", nrow(gcimp.samples),"\n",
              "GCIMP Low samples ", nrow(gcimp.low.samples),"\n",
              "GCIMP High samples ", nrow(gcimp.high.samples),"\n",
              "-----------------------------------------------")
message(msg)
#-----------------------------------
# 2 - Get data
# ----------------------------------
#-----------------------------------
# 2.1 - DNA Methylation
# ----------------------------------
for(cancer in c("LGG","GBM")){
  query <- GDCquery(project = paste0("TCGA-",toupper(cancer)),
                    data.category = "DNA Methylation",
                    sample.type = "Primary solid Tumor",
                    platform = "Illumina Human Methylation 450",
                    barcode =  gcimp.samples$patient)
  tryCatch({
    GDCdownload(query, chunks.per.download = 20)
  }, error = function(e) {
    GDCdownload(query, chunks.per.download = 5)
  })
  assign(tolower(cancer),GDCprepare(query, save = FALSE))
}
met <- SummarizedExperiment::cbind(lgg, gbm)
save(met, file = "gcimp_met.rda")
rm(lgg,gbm)
#-----------------------------------
# 2 - Expression
# ----------------------------------
for(cancer in c("LGG","GBM")){
  query <- GDCquery(project = paste0("TCGA-",toupper(cancer)),
                    data.category = "Transcriptome Profiling",
                    data.type = "Gene Expression Quantification",
                    workflow.type = "HTSeq - FPKM-UQ",
                    sample.type = "Primary solid Tumor",
                    barcode =  gcimp.samples$patient)
  tryCatch({
    GDCdownload(query, chunks.per.download = 20)
  }, error = function(e) {
    GDCdownload(query, chunks.per.download = 5)
  })
  assign(tolower(cancer), GDCprepare(query, save = FALSE))
}
exp <- SummarizedExperiment::cbind(lgg, gbm)
save(exp, file = "gcimp_exp.rda")
rm(lgg,gbm)

genome <- "hg38"
file <- "hg38_gcimp_mae.rda"
distal.probe <- get.feature.probe(genome = genome, met.platform = "450K",feature = NULL)
if(file.exists(file)) {
  mae <- get(load(file))
} else {
  mae <- createMAE(exp = "gcimp_exp.rda",
                   met = "gcimp_met.rda",
                   genome = genome,
                   met.platform = "450K",
                   linearize.exp = TRUE,   # Log2(exp + 1)
                   filter.probes = distal.probe,
                   TCGA = TRUE,
                   save = TRUE,
                   save.filename = file)
}
cores <- parallel::detectCores()/2
direction <- c("hyper","hypo")
group1 <- "G-CIMP-high"
group2 <- "G-CIMP-low"
group.col <- "subtype_Supervised.DNA.Methylation.Cluster"

for (j in direction){
  print(j)
  dir.out <- paste0("gcimp_high_vs_low_",genome,"/",j)
  dir.create(dir.out, recursive = TRUE)
  #--------------------------------------
  # STEP 3: Analysis                     |
  #--------------------------------------
  # Step 3.1: Get diff methylated probes |
  #--------------------------------------
  Sig.probes <- get.diff.meth(data       = mae,
                              group.col  = group.col,
                              group1     = group1,
                              group2     = group2,
                              minSubgroupFrac = 1,
                              cores      = cores,
                              dir.out    = dir.out,
                              diff.dir   = j,
                              pvalue     = 0.01)

  #-------------------------------------------------------------
  # Step 3.2: Identify significant probe-gene pairs            |
  #-------------------------------------------------------------
  # Collect nearby 20 genes for Sig.probes
  nearGenes <- GetNearGenes(data  = mae,
                            probe = Sig.probes$probe,
                            cores = cores)

  pair <- get.pair(data       = mae,
                   nearGenes  = nearGenes,
                   group.col  = group.col,
                   group1     = group1,
                   group2     = group2,
                   permu.dir  = paste0(dir.out,"/permu"),
                   dir.out    = dir.out,
                   diff.dir   = j,
                   minSubgroupFrac = 1,
                   cores      = cores,
                   mode       = "supervised",
                   label      = j,
                   permu.size = 10000,
                   raw.pvalue = 0.001)
  Sig.probes.paired <- unique(read.csv(paste0(dir.out,
                                              "/getPair.",j,
                                              ".pairs.significant.csv"),
                                       stringsAsFactors=FALSE)[,1])


  #-------------------------------------------------------------
  # Step 3.3: Motif enrichment analysis on the selected probes |
  #-------------------------------------------------------------
  if(length(Sig.probes.paired) > 10 ){
    #-------------------------------------------------------------
    # Step 3.3: Motif enrichment analysis on the selected probes |
    #-------------------------------------------------------------
    enriched.motif <- get.enriched.motif(probes  = Sig.probes.paired,
                                         dir.out = dir.out,
                                         data    = mae,
                                         label   = j)
    motif.enrichment <- read.csv(paste0(dir.out,
                                        "/getMotif.",j,
                                        ".motif.enrichment.csv"),
                                 stringsAsFactors=FALSE)
    if(length(enriched.motif) > 0){
      #-------------------------------------------------------------
      # Step 3.4: Identifying regulatory TFs                        |
      #-------------------------------------------------------------
      print("get.TFs")

      TF <- get.TFs(data           = mae,
                    enriched.motif = enriched.motif,
                    group.col      = group.col,
                    group1         = group1,
                    group2         = group2,
                    dir.out        = dir.out,
                    minSubgroupFrac = 1,
                    diff.dir       = j,
                    cores          = cores,
                    mode           = "supervised",
                    label          = j)
      TF.meth.cor <- get(load(paste0(dir.out,
                                     "/getTF.",j,
                                     ".TFs.with.motif.pvalue.rda")))
      save(mae,TF, enriched.motif, Sig.probes.paired,
           pair, nearGenes, Sig.probes, motif.enrichment,
           TF.meth.cor,
           file = paste0(dir.out,"/ELMER_results_",j,".rda"))
    }
  }
}
for (j in direction){
  print(j)
  dir.out <- paste0("gcimp_high_vs_low_",genome,"/",j)
  pairs <- readr::read_csv(file.path(file.path(dir.out,paste0("getPair.",j,".pairs.significant.csv"))))
  message("Scatter plots")
  scatter.plot(mae,
               byPair = list(probe = pairs$Probe, gene = pairs$GeneID),
               category = "subtype_Supervised.DNA.Methylation.Cluster",
               save = TRUE,
               lm_line = TRUE,
               dir.out = paste0(dir.out,"/significant_pairs_plot"))
  message("Schematic plot for probes")
  schematic.plot(data = mae,
                 group.col = "subtype_Supervised.DNA.Methylation.Cluster",
                 pair = pairs,
                 dir.out = paste0(dir.out,"/schematic_plots"),
                 save = TRUE,
                 byProbe = unique(pairs$Probe))
  message("Schematic plot for genes")
  # by gene ID and adding annotation tracks from StateHub               
  schematic.plot(data = mae,
                 group.col = "subtype_Supervised.DNA.Methylation.Cluster",
                 pair = pairs,
                 dir.out = dir.out,
                 save = TRUE,
                 byGeneID = unique(pairs$GeneID))
}

4 Session Information

## R version 3.4.3 (2017-11-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
## LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0
##
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
##  [1] parallel  stats4    grid      stats     graphics  grDevices utils
##  [8] datasets  methods   base
##
## other attached packages:
##  [1] bindrcpp_0.2               GenomicRanges_1.30.1
##  [3] GenomeInfoDb_1.14.0        IRanges_2.12.0
##  [5] S4Vectors_0.16.0           BiocGenerics_0.24.0
##  [7] funciVar_0.1.4             ggplot2_2.2.1.9000
##  [9] ComplexHeatmap_1.17.1      knitr_1.17.20
## [11] MultiAssayExperiment_1.4.4 dplyr_0.7.4
## [13] ELMER_2.3.6                ELMER.data_2.3.2
## [15] DT_0.2
##
## loaded via a namespace (and not attached):
##   [1] shinydashboard_0.6.1          R.utils_2.6.0
##   [3] RSQLite_2.0                   AnnotationDbi_1.40.0
##   [5] htmlwidgets_0.9               BiocParallel_1.12.0
##   [7] devtools_1.13.4               DESeq_1.30.0
##   [9] munsell_0.4.3                 codetools_0.2-15
##  [11] withr_2.1.1.9000              colorspace_1.3-2
##  [13] BiocInstaller_1.28.0          Biobase_2.38.0
##  [15] highr_0.6                     rstudioapi_0.7
##  [17] labeling_0.3                  GenomeInfoDbData_1.0.0
##  [19] mnormt_1.5-5                  hwriter_1.3.2
##  [21] KMsurv_0.1-5                  bit64_0.9-7
##  [23] rprojroot_1.3-1               downloader_0.4
##  [25] biovizBase_1.26.0             ggthemes_3.4.0
##  [27] EDASeq_2.12.0                 R6_2.2.2
##  [29] doParallel_1.0.11             locfit_1.5-9.1
##  [31] AnnotationFilter_1.2.0        bitops_1.0-6
##  [33] reshape_0.8.7                 DelayedArray_0.4.1
##  [35] assertthat_0.2.0              scales_0.5.0.9000
##  [37] nnet_7.3-12                   gtable_0.2.0
##  [39] ensembldb_2.2.0               rlang_0.1.6
##  [41] genefilter_1.60.0             cmprsk_2.2-7
##  [43] GlobalOptions_0.0.12          splines_3.4.3
##  [45] rtracklayer_1.38.2            lazyeval_0.2.1
##  [47] acepack_1.4.1                 dichromat_2.0-0
##  [49] selectr_0.3-1                 broom_0.4.3
##  [51] checkmate_1.8.5               yaml_2.1.16
##  [53] reshape2_1.4.3                snpStats_1.28.0
##  [55] GenomicFeatures_1.30.0        backports_1.1.2
##  [57] httpuv_1.3.5                  Hmisc_4.0-3
##  [59] RMySQL_0.10.13                tools_3.4.3
##  [61] psych_1.7.8                   RColorBrewer_1.1-2
##  [63] Rcpp_0.12.14                  plyr_1.8.4
##  [65] base64enc_0.1-3               progress_1.1.2
##  [67] zlibbioc_1.24.0               purrr_0.2.4
##  [69] RCurl_1.95-4.8                prettyunits_1.0.2
##  [71] ggpubr_0.1.6                  rpart_4.1-11
##  [73] pbapply_1.3-3                 GetoptLong_0.1.6
##  [75] zoo_1.8-0                     SummarizedExperiment_1.8.1
##  [77] ggrepel_0.7.0                 cluster_2.0.6
##  [79] magrittr_1.5                  data.table_1.10.4-3
##  [81] circlize_0.4.3                survminer_0.4.1
##  [83] ProtGenerics_1.10.0           matrixStats_0.52.2
##  [85] aroma.light_3.8.0             hms_0.4.0
##  [87] mime_0.5                      evaluate_0.10.1
##  [89] xtable_1.8-2                  XML_3.98-1.9
##  [91] gridExtra_2.3                 shape_1.4.3
##  [93] testthat_2.0.0                compiler_3.4.3
##  [95] biomaRt_2.34.1                tibble_1.3.4
##  [97] R.oo_1.21.0                   htmltools_0.3.6
##  [99] Formula_1.2-2                 tidyr_0.7.2
## [101] geneplotter_1.56.0            DBI_0.7
## [103] matlab_1.0.2                  ShortRead_1.36.0
## [105] Matrix_1.2-12                 readr_1.1.1
## [107] R.methodsS3_1.7.1             Gviz_1.22.2
## [109] bindr_0.1                     pkgconfig_2.0.1
## [111] km.ci_0.5-2                   GenomicAlignments_1.14.1
## [113] foreign_0.8-69                plotly_4.7.1
## [115] xml2_1.1.1                    roxygen2_6.0.1
## [117] foreach_1.4.4                 annotate_1.56.1
## [119] XVector_0.18.0                rvest_0.3.2
## [121] stringr_1.2.0                 VariantAnnotation_1.24.2
## [123] digest_0.6.13                 ConsensusClusterPlus_1.42.0
## [125] Biostrings_2.46.0             rmarkdown_1.8
## [127] TCGAbiolinks_2.7.4            survMisc_0.5.4
## [129] htmlTable_1.11.0              edgeR_3.20.2
## [131] curl_3.1                      shiny_1.0.5
## [133] Rsamtools_1.30.0              commonmark_1.4
## [135] rjson_0.2.15                  nlme_3.1-131
## [137] jsonlite_1.5                  viridisLite_0.2.0
## [139] limma_3.34.4                  BSgenome_1.46.0
## [141] lattice_0.20-35               httr_1.3.1
## [143] survival_2.41-3               interactiveDisplayBase_1.16.0
## [145] glue_1.2.0                    iterators_1.0.9
## [147] bit_1.1-12                    stringi_1.1.6
## [149] blob_1.1.0                    AnnotationHub_2.10.1
## [151] latticeExtra_0.6-28           memoise_1.1.0

ELMER - GCIMP analysis - G-CIMP-high vs G-CIMP-low

Tiago Chedraoui Silva, Thais Sabedot Serraf, Ben P. Berman, Houtan Noushmehr

2017-12-25

1 Summary

1.1 Groups information

1.2 Samples

1.3 Arguments

1.4 Summary results

1.4.1 Probes hypermethylated in G-CIMP-high vs G-CIMP-low

2 Glioma analysis: G-CIMP-high vs G-CIMP-low

2.1 Probes hypermethylated in G-CIMP-high vs G-CIMP-low

2.1.1 Significant anti-correlated pairs of gene-probes

2.1.1.1 Statehub: Chromatin state evaluation

2.1.2 Motif enrichment analysis

2.1.3 TF analysis

2.1.3.1 Top potential TF

2.1.3.2 potential TF family

2.1.3.3 potential TF subfamily

2.1.3.4 TF plots

2.2 Probes hypomethylated in G-CIMP-high vs G-CIMP-low

3 Code

4 Session Information

ELMER - GCIMP analysis - G-CIMP-high vs G-CIMP-low

Tiago Chedraoui Silva, Thais Sabedot Serraf, Ben P. Berman, Houtan Noushmehr

2017-12-25

1 Summary

1.1 Groups information

1.2 Samples

1.3 Arguments

1.4 Summary results

1.4.1 Probes hypermethylated in G-CIMP-high vs G-CIMP-low

2 Glioma analysis: G-CIMP-high vs G-CIMP-low

2.1 Probes hypermethylated in G-CIMP-high vs G-CIMP-low

2.1.1 Significant anti-correlated pairs of gene-probes

2.1.1.1 Statehub: Chromatin state evaluation

2.1.2 Motif enrichment analysis

2.1.3 TF analysis

2.1.3.1 Top potential TF

2.1.3.2 potential TF family

2.1.3.3 potential TF subfamily

2.1.3.4 TF plots

2.1.3.4.1 TF plots

2.1.3.4.1.1 LHX2

2.1.3.4.1.2 PO6F2

2.1.3.4.1.3 RX

2.1.3.4.1.4 JUNB

2.1.3.4.1.5 GBX2

2.1.3.4.1.6 NOTO

2.1.3.4.1.7 DLX3

2.1.3.4.1.8 SOX10

2.1.3.4.1.9 MIXL1

2.1.3.4.1.10 SIX2

2.1.3.4.1.11 LBX2

2.1.3.4.1.12 BATF

2.1.3.4.1.13 EVX1

2.1.3.4.1.14 ZBT49

2.1.3.4.1.15 SIX1

2.1.3.4.1.16 SOX4

2.1.3.4.1.17 DLX2

2.1.3.4.1.18 FOXQ1

2.1.3.4.1.19 SOX9

2.1.3.4.1.20 HXB3

2.1.3.4.1.21 DLX5

2.1.3.4.1.22 HXD3

2.1.3.4.1.23 LHX9

2.1.3.4.1.24 VAX2

2.1.3.4.1.25 EMX2

2.1.3.4.1.26 PAX4

2.1.3.4.1.27 GBX1

2.1.3.4.1.28 FOSB

2.1.3.4.1.29 JUND

2.1.3.4.1.30 JUN

2.1.3.4.1.31 SOX13

2.1.3.4.1.32 TWST1

2.1.3.4.1.33 HME2

2.1.3.4.1.34 FOSL2

2.1.3.4.1.35 ELF3

2.1.3.4.1.36 NKX61

2.1.3.4.1.37 SOX2

2.1.3.4.1.38 NFIA

2.1.3.4.1.39 SOX10

2.1.3.4.1.40 BARX1

2.1.3.4.1.41 EMX1

2.1.3.4.1.42 HXA2

2.1.3.4.1.43 DLX4

2.1.3.4.1.44 NF2L2

2.1.3.4.1.45 BACH2

2.1.3.4.1.46 PO4F2

2.1.3.4.1.47 ATF2

2.1.3.4.1.48 TEAD2

2.1.3.4.1.49 FOXC2

2.1.3.4.1.50 VSX2

2.1.3.4.1.51 CEBPA

2.1.3.4.1.52 SOX9

2.1.3.4.1.53 EHF

2.1.3.4.1.54 FOS

2.1.3.4.1.55 SOX3

2.1.3.4.1.56 NFE2

2.1.3.4.1.57 NKX62

2.1.3.4.1.58 HXA13

2.1.3.4.1.59 SOX18

2.1.3.4.1.60 OLIG2